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Abstract 

This work combines three paradigms of image processing: i) the total 
variation approach to denoising, ii) the superior structure of hexagonal lat- 
tices, and iii) fast and exact graph cut optimization techniques. Although 
isotropic in theory, numerical implementations of the BV seminorm invari- 
ably show anisotropic behaviour. Discretization of the image domain into a 
hexagonal grid seems perfectly suitable to mitigate this undesirable effect. To 
this end, we recast the continuous problem as a finite-dimensional one on an 
arbitrary lattice, before focussing on the comparison of Cartesian and hexag- 
onal structures. Minimization is performed with well-established graph cut 
algorithms, which are easily adapted to new spatial discretizations. Apart 
from producing minimizers that are closer in the i-^ sense to the clean image 
for sufficiently high degrees of regularization, our experiments suggest that 
the hexagonal lattice also allows for a more effective reduction of two major 
drawbacks of existing techniques: metrication artefacts and staircasing. For 
the sake of practical relevance we address the difficulties that naturally arise 
when dealing with non-standard images. 

Keywords, image denoising, variational methods, total variation regularization, 
hexagonal lattice, graph cuts 



1 Introduction 

1.1 Total Variation Denoising 

Denoising is a prominent image processing task, and total variation (TV) regular- 
ization is one particular approach. It consists in recovering a noise-free image from 
perturbed data / G by finding the solution of an optimization problem, 

which can be written in general form as 

min A||m-/||2„+TV(m). {TVL") 

u^BV {il) 

We restrict our attention to bounded domains Q C and only consider scalar 
functions f,u > 0, which have an immediate interpretation as common grey scale 
images. Possible minimizers of ( TVL^| are required to be of bounded variation. 



in other words, the space BV{Q) ^ {u e L\Q) : TV(i^) < oo} is adopted for 
the modelling of images. For Sobolev functions u G VK^'^(n) C BV{Q) the total 
variation equals \ Vu\. However, it is usually defined via duality as 

TY{u) = sup|y i^div^dx : g G Co°°(Q,R^), \g{x)\ < 1, x G (1.1) 

to include discontinuous functions, which is actually one of the main virtues of this 
denoising technique. Equipped with the norm |H|^^ = IMIl1(Q) +TV(-) the space 
of functions of bounded variation is Banach. We refer to [4j [28] for theoretical 



background on BV(Q), to [T4l[T2l[5T] for more details on its role in image analysis 
and to 37| for a concise comparison to other image filters. 



Returning to (TVL""), we only consider the two popular cases a G {1,2}, 
but other choices are possible. For a = 2 the minimization problem becomes the 
unconstrained Rudin-Osher-Fatemi model [49l[TJ[T6], which was originally designed 
to remove Gaussian noise. An data term is also widely used and seems to be 
more suitable for impulsive noise [3] 1421 1431 I18 L From an analytical point of 
view, one major difference between the two models is that the former is strictly 
convex, thus having a unique minimizer in BV(Q), while the latter is only convex, 
allowing for multiple minimizers. The weighting parameter A > controls the 
trade-off between data fidelity and degree of regularization. Due to its ability to 
preserve sharp edges, the TV approach to image analysis is considered a substantial 
improvement over previous least squares methods using quadratic regularization 
terms. However, it still suffers from several shortcomings, such as staircasing and 
loss of texture [39] . 

Over the last two decades a great number of different numerical methods for 



solving (TVL"") have been developed. In the words of [17], they can be divided 
into two classes: (i) those that "optimize" before discretizing and (ii) those that do 
it in reverse order. A common approach that belongs to the first group consists 
in solving an artificial time evolution of the associated Euler- Lagrange equation 
with some difference scheme. The second class deals with a discrete version of the 
original functional, thus solves a finite dimensional minimization problem. Graph 
cuts, which we chose to employ in this work, belong to the latter. 

1.2 Graph Cuts 

Essentially, graph cut minimization techniques make use of two equivalences, both 
of which were identified several decades ago: first, the one between certain dis- 
crete functions and cut functions of s-t-graphs [3^ |46j and, second, that between 
minimum cuts and maximum flows on such graphs [27]. In order to minimize the 
objective function, the edge weights in the graph are chosen so that there is a 
one-to-one correspondence between costs of cuts and values of the function. Since 
maximum flows can be found in low-order polynomial time with small constants, 
solutions can be computed very efficiently |2][8]. Digital images naturally fit into 
the graph cuts framework, which is why the computer vision community have ex- 
ploited their benefits for a wide variety of problems. For a brief introduction to 
graph cuts see [9]. For an in-depth treatment with numerous extensions, applica- 
tions and a comprehensive bibliography we refer to the recent handbook ^6 . 

Concerning the problem of TV denoising, there are several works that employed 



graph cuts to compute exact minima of discrete versions of (TVL^ 



>'^\Up- /pT + ^^pq\Up - Uq\ (1.2) 

p p,q 

in similar ways — see e.g. i33l 1131 l22l 1151 129j (in chronological order). Variables p 
and q denote sites of a discretization of domain Q and ojpq are weights according 



to some discrete TV (cf. Sec. 2.3). Using a level set formulation, energy (1.2) is 
decomposed into a sum of binary functions, which in turn are minimized separately. 
The fact that, on the one hand, TV admits a complete decoupling of levels of u and 
that, on the other hand, flow information from one level can be reused for another, 



makes solving (TVL"" ) highly efficient. Minimization of the multi-label energy ( 1.2 ) 
in one large sweep with the method from [34] was reported to be considerably less 
efficient, in terms of both running time and memory consumption 13 . 

Although the TV, as defined in ( |1.1[ ), is isotropic, that is it assigns the same 
values to rotated versions of an image, it is well-known that numerical implementa- 
tions tend to be anisotropic [111 [T2] . Due to the two-directional nature of the pixel 
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grid (and thus of the flow network), graph cut approaches are no exception. One 
way of mitigating this undesirable effect is to consider grids with larger neighbour- 
hood structures. Another possibility is to abandon the square lattice altogether, 
in favour of a hexagonal one. 



1.3 Hexagonal Grids 

That a rectangular lattice might not be the best choice for the sampling and pro- 
cessing of two-dimensional signals has been recognized at least half a century ago 
|45l|38]. Since then, researchers who investigated the use of hexagonally arranged 
pixel configurations almost unanimously conclude that they should be preferred 
over their more common counterparts for a wide variety of applications, such as 
edge detection, shape extraction or image reconstruction (see e.g. the textbook 
[40] . survey articles [H [31] or, in chronological order, [Ml Ell [H [25]). The 
virtues of a hexagonal grid — dense packing, better connectivity, higher angular 
resolution — have also been acknowledged in the field of morphological image anal- 
ysis [551 152j as well as for novel imaging frameworks 2A^ . They can be traced back 
to a few closely related geometric optimality properties [201 [30] . Inspired by these 
and their frequent occurrence in nature (insect and human visual systems, honey- 
combs) hexagonal arrays are heavily used in the area of optical technology (e.g. for 
large telescope mirrors). Finally, let us remark that hexagonal pixel arrangements 
are already implemented in robotics vision systems}^ largescale media display^ 
and electronic paper 32 . 

The obvious drawback of a hexagonal grid is its non-alignment with Cartesian 
coordinates, which means that, in contrast to a square lattice, there is no straight- 
forward a) calculus, b) extension to higher dimensions and c) implement ability. 
This, at least to some extent, explains the present lack of hardware devices for 
efficiently capturing, storing and visualizing hexagonally sampled images, which 
must be overcome by means of software (see ^40^ for a recently proposed hexagonal 
image processing framework). In Figs, [l] and |2] hexagonally sampled images are 
opposed to common square pixel images. Due to the sufficiently high resolution, 
the difference between the two sampling schemes seems only marginal in Fig. [l] 
In Fig. [2] which features smooth intensity changes, the better visual effects of 
a hexagonal pixel configuration becomes more apparent. Throughout this note, 
hexagonal pixels are visualized using so-called hyperpixels, i.e. groupings of square 
pixels, that are approximately hexagonal [40] . 



1.4 Outline 



When continuous problems such as (TVL"" ) are brought to finite dimensions, the 



domain Q is usually transformed into a discrete set of points arranged in an or- 
thogonal way. Consequently, functions defined on this lattice inevitably inherit its 
strongly anisotropic character. However, apart from an easy calculus, there is ac- 
tually no good reason for choosing a rectangular grid, in particular when isotropy 
is an important property of the underlying problem as is the case for image pro- 
cessing. The fact that a hexagonal grid features three natural directions instead of 
only two can be expected to mitigate undesirable discretization errors. 

To be able to exploit the advantages of the hexagonal grid, we first formulate 



a finite-dimensional version of (TVL" ) on an arbitrary lattice. The corresponding 
discrete functional is formally derived in Section[2] In Sec. |2.1| the discretization of 
domain Q is discussed, that is, some background on lattices and neighbourhoods is 



provided. From that, the definition of discrete images is straightforward (Sec. 2.2 ) 



Sec. 12.31 is devoted to the discretization of the TV functional. In Sec. 12.41 we state 



""^http:/ /centeye.com/ 

^http:/ /www. smartslab.co.uk/ 
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Figure 1: Contrast-enhanced version of the original Shepp-Logan head phantom |53| , sampled 
to a 100 X 100 square pixel grid (right), and to a hexagonal one with the same resolution (left). 





Figure 2: Radially symmetric cosine, cos ( ^ 45^0 ^ i^jU) ^ [0? 100] X [—100,0], scaled to the 
interval [0, 255] and sampled to both a 54 x 54 square pixel grid (right) and a hexagonal grid 
of the same sampling density (left). 



our central problem, an "optimal" finite-dimensional version of problem (TVL^^ 



on an arbitrary lattice, and in Sec. |2.5| we briefly examine the relations between 
the continuous and discrete functionals. Finally, in Sec. [3] we present experimental 
results for hexagonal pixel grid discretizations compared to traditional square ones. 
We conclude this paper with a brief discussion of our results in Sec. [4] 

We emphasize that the main contribution of the present work is neither a novel 
image restoration algorithm nor a theoretical result, but rather a qualitative and 
purely empirical investigation of different two-dimensional discretization schemes 
for the problem of TV denoising, focussing on the comparison of square pixel 
structures with hexagonal ones. 



Figure 3: Horizontally aligned regular hexagonal lattice with Voronoi cell tiling. 



2 The Discrete Setting 
2.1 Planar Lattices 

A planar lattice is a regular discrete subset of the Euclidean plane. Any basis 
B = {bi, 1)2) of generates a lattice A(B) in the following way 

K{B) = [Bp -.pel?]. (2.1) 

Two bases Bi,B2 define the same lattice A{Bi) = A(B2), if and only if one can 
be obtained from the other with elementary integer column operations, i.e. by 
swapping columns, multiplication of columns by —1 and adding integer multiples of 
one column to another. In other words, there exists a unimodular matrix U G R^^^, 
an integer matrix with integer inverse, so that Bi = B2U. 

Example 1. Let us give two examples that will be used throughout this paper. 
First, if B = / is the identity matrix, A{B) = is the usual unit square lattice. 
Second, if the basis vectors 61, 62 have the same length and draw an angle of | or 
then A(B) becomes a regular hexagonal lattice. For the sake of simplicity we 
only consider horizontally aligned hexagonal lattices, and set accordingly 

Fig. [3] shows a horizontally aligned hexagonal lattice. A 
2.1.1 Lattice Density and Tilings 

From now on, we restrict ourselves to bases B consisting of unit vectors, and 
introduce an additional parameter d > controlling the distance between points 
of lattice A(dB). A discretization of the image domain Q is then defined as 

:=nnA{dB). (2.3) 

We set N = and, in order to simplify notations, we will identify a lattice 

point dBp G ^% simply by its "index" p. In general, d does not equal the min- 
imal distance between two lattice points, although it does in the two interesting 
cases A{dl) and A{dH). Parameter d does, however, equal the side length of the 
lattice's fundamental parallelogram U{dB) = {dBx : X G [0, 1)^}. The parallelo- 
gram's area |det(ii?| is an important quantity, since it is inversely proportional 
to lattice density: there is one lattice point per parallelogram. From it we can 
deduce the phenomenon, which is sometimes referred to as the "denser packing" of 
hexagonal lattices: 

/o 

\det dH\ = ^d^ <d^ = \det dl\. (2.4) 
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This implies that A(dH) has a higher density than A(dl) and that, in general, 
> win hold. Obviously, if we want to compare different lattices, it is 



important that they are equally dense. Therefore, due to (2.4), we have to shrink 
the square lattice by a factor of ^/2, or stretch the hexagonal one by y^^/Vs, in 
order to get a square lattice with the same density as a hexagonal one — and, thus, 
a discretization of Q with approximately the same amount of lattice points. 

Associated to every lattice, there is one monohedral tiling, i.e. a partition of the 
plane consisting of congruent shapes, which will be of use later on. It is made up 
of Voronoi cells (cf. Fig.jsj). The Voronoi cell associated to one point p belonging 
to some discrete set A C is defined as the set of all points that are closer to p 
than to any other point from A 

V{p) = {xeM.^ :\x-p\< \x-q\,\/qeA\{p}}. (2.5) 

The area of a Voronoi cell associated to lattice A{B) is equal to that of the funda- 
mental parallelogram |det B\, since Voronoi cells and translations of the fundamen- 
tal parallelogram both form monohedral tilings of the plane with one lattice point 
per tile. In the following we will mainly be interested with that part of a Voronoi 
cell that lies in Q, i.e. with Vp := V{p)nQ. For more details on (hexagonal) lattices 
we refer to I20ll40l. 



2.1.2 Neighbourhoods 



There are several equivalent ways of defining a neighbourhood relation between 
lattice points. We define it (as in [7j) by a set of vectors pointing in the direction 
of neighbouring sites 

A/-^ = {^, GR^:i = l,...,z/}. (2.6) 

These vectors are the "directions" of the lattice and, since neighbourship is a sym- 
metric relation, by Vi we actually mean ±Vi. It is also clear that these vectors are 
linear combinations of the lattice's basis vectors. Two sites p and q are neighbours 
according to A/'b, if dB{p — q) equals some Vi G Mb- We will indicate this with the 
shorthand p ^ q, when there is no confusion about the neighbourhood system. 



For the square lattice at least three neighbourhood systems (Fig. 4a) have been 
proposed and already used in the context of TV minimization |29j : 



N't = {ei,e2}, 

jVf = jVf U {ei + 62, 62 - 6i} , 
A/-/6 = A/f U {261 + 62, 61 + 262, 262 ' 



(2.7) 



61, 62 - 26i} . 



Here, Ci are the canonical basis vectors and we set d — 1 for simplicity. For the 



hexagonal lattice two natural neighbourhoods are (Fig. 4b): 
J^H - {hi,h2, h2 - hi} , 



^A/'^U{/ii + /i2,2/i2 -hi,h2 -2hi}. 



(2.^ 



The important point here is that \h2 — hi\ = \hi\ — \h2\. A regular hexagonal 
lattice is the only planar lattice that permits each of its points to have six nearest 
neighbours instead of only four or even less. 



2.2 Generalized Discrete Images and their Levels 

Throughout this note we regard a generalized discrete image u to be a quantized 
and sampled version of G BV{Q). That is, u is rounded to a discrete set of 
grey values £ C M, and evaluated at sites of lattice Q%. According to the notation 
introduced in the previous section, we write 

u {up)p^n% ^ Up := [u{p)\, (2.9) 
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(a) Neighbours of the central black 
point according to Aff (thick lines), 
(thick and fine) and J\f]^ 
(thick, fine and dashed). 
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(b) Neighbours of the central black point 
according to M% (thick lines), Af}^ 
(thick and fine) 



Figure 4: Square lattice (left) and hexa gona l lattice (right) with neighbourhoods from \2.1\ 

and (|2.8|), respectively. 



Table 1: Super level set characteristics of a discrete image. 
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where [-J stands for quantization |^ Such an image is "general" in the sense that 
for a rectangular domain Q and B — I \i reduces to the widely used matrix 
representation of digital images. Additionally, in order to simplify notations, C is 
assumed to be a set of L consecutive integers starting from zero, i.e. C — {0, . . . , L — 
1}, where and L — 1 are visualized as black and white, respectively, and values 
in between are represented by different shades of grey (cf. [5]). 
For / G M the super level sets of u are defined as 

Ei^{x^VL'. u{x) > I}. (2.10) 

We will denote the characteristic function of a super level set Ei of u simply 
by x^ x\^) — I5 if u(x) > I and x^(^) — otherwise. Image levels are 
analogously defined in the discrete setting. To each discrete image u, L binary 
images x' — (Xp)pGO^ 5 ^ ^ ^1 ^ire associated, where Xp equals 1, ii Up is above level 
/, and if Up is below or equal to /. Note that there are actually only L — 1 relevant 
levels, since x^~^ = contains no information. This family of binary images is 
monotone decreasing, i.e. it satisfies 

l>k ^ Xp<Xp, P^^i- (2.11) 

In other words, for each lattice point p they form a monotone decreasing sequence 
as depicted in Tab. [l] Knowing all levels of an image is equivalent to knowing the 
image itself. Thus u can be reconstructed from its levels by Up — Xp — min{/ : 
Xp = 0}. 

2.3 Discrete Total Variation 

In this section we demonstrate how to discretize the TV functional in a reasonable 
way, by combining the coarea and Cauchy-Crofton formulae. While the former 

^Since u G BV{Q) C L^{Q), we think of u{p) as lim^^Q_(_ (r^Tr)"""^ n, where Br{p) is 

the disk of radius r centred at p. 
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states that the TV can be computed by summing level contours, the latter provides 
a sound way of estimating those contours from function values on a lattice. This is 
the common approach for graph cut approaches to TV denoising [22] [29] . Below, 
we demonstrate that this approach is easily generalized to arbitrary lattices and, 
by doing so, make a further point in favour of the hexagonal grid. 
For u G BV{Q) the coarea formula reads 



/_ 



TV(m)= / TV{x)dl, (2.12) 



where is the characteristic function of the set Ei as defined in (2.10) [31 128]. 
The quantity PeT{Ei) := TV(x^) is finite for almost all / G M and is commonly 
called perimeter of the level set. The coarea formula thus states that the TV of 
a BF-function can be computed by summing the perimeters of all its super level 
sets|^ In the following we denote the n-dimensional Hausdorff measure of sets by 
vertical bars. From the context it should be clear whether the number of elements 
(n = 0), length (n = 1) or area (n = 2) is meant. 

2.3.1 The Cauchy-Crofton Formula 

The Cauchy-Crofton formula is an elegant result from integral geometry, relating 
the length of a plane curve to the measure of straight lines intersecting it [23] . 
Inspired by the work of |7| , we utilize it to estimate from u the level set perimeters 
of 2/ in a sound way. 

Consider a straight line /c in . By identifying it with its signed distance from 
the origin p and the angle which it draws with the x-axis, line k can be regarded 
as a point in the set AC = R x [0, tt)]^ Let C be a rectifiable plane curve C and nc 
the function mapping each line to the number of times it intersects with C. The 
Cauchy-Crofton formula states that 

\C\ = l [ nc{p,<P)dpd<p. (2.13) 



One of the virtues of equation (2.13) is that it can also be used to approximate 
the length of C: by replacing the double integral on the right-hand side with a 
partial sum, which corresponds to equipping the plane with count ably many lines, 
a reasonable estimate for the curve length can be obtained. More precisely, we have 
m families of lines Fi, each characterized by its angle and inter-line distance 
Api, i.e. 

F,(0,,ApO = {(jAp,,0O :jeZ}, i = l,...,m. (2.14) 
In addition, we assume = 01 < 02 < • • • < 0m < vr and set A0i = 0i+i — 0^ for 



i < m — 1, and A0 771 — TT (pm • Then, from (^2.13fc we get 



1^1 - 7;^^ric{jAp^,(l)^)/\p^/\(l)^^^nc{jAp^,(l)^)uj^. (2.15) 



=1 jez 



Here, we used the shorthand uji = ApiA0i/2. This approximation formula tells 
us for each family Fi to count its intersections with C and to multiply the result 



with half the area of the rectangle in /C corresponding to Fi. Moreover, (2.15) is 



^In the case where Ei has a boundary of class (up to a set which is negligible with respect 
to the one-dimensional Hausdorff measure) one has \dEi \ = Per(£^;), that is, the perimeter of Ei 
indeed equals the length of its boundary. In general, however, for u G BV{yi) the topological 
boundary is not the appropriate object of study and must be replaced with the (measure- 
theoretic) reduced boundary Ei, a dense subset of dEi [4] I28|. 

^To make the assignment unique, we set (p, ^) to be the line x = p. Note that the Cauchy- 
Crofton formula is usually stated using normal parameters (p,(f), i.e. its unsigned distance p 
to the origin and the angle ip which is drawn by the a:-axis and a line normal to k |23| . 



8 



a reasonable approximation, since it converges to the true length of C as m ^ oo 
and sup - Ap^ 0, sup - A(/)i [23] . 

Consider a grid graph embedded in the same plane where C lies, consisting of a 
set of nodes A{dB) and edges according to some neighbourhood A/"^. Such a graph 
naturally induces m — ^ families of straight lines Fi. Furthermore, we assign 
weights (jJi from ( 2.15| ) to edges, depending on the family to which they belong. 



Now, by summing up the weights of all those edges that have been intersected by 



C, we get exactly the right hand side of (2.15), provided the lattice is sufficiently 



fine so that each edge is intersected by C at most once. Due to its higher angular 
resolution, a hexagonal grid graph can be expected to yield a better approximation 
than a square one, since it naturally equips the plane with more families of lines. 

Example 2 (Weights for A/"^, M}^). In the case of a regular hexagonal lattice 
A(dH) with neighbour distance d > and neighbourhood A/"^, we have m = 3 
families of lines at angles 0^ = {i — 1)^ and line distance Api — Ap = dVS/2. 
Therefore, we get cj^ = a; = VSf^d. 



For JV}^ (cf. Fig. 4b) we have m = 6 families at angles = (i — and line 
distances 



Ap^ 



2 ^ 



for i G {1,3,5} 
for i G {2,4,6}. 



Accordingly, we get two different weights: one for edges between nearest neighbours 
standing at distance and one for edges connecting sites at distance 



^d, for {1,3,5} 
^d, for {2,4,6}. 



. 24"- 

The weights for Nf , Nf and N}^ can be found in [29]. A 
2.3.2 TV Estimation on Arbitrary Lattices 

Consider a continuous signal u G BViO)^ which is known to us onl y at sampling 



points Q^B^ as a discrete image u on an arbitrary lattice (cf. Sec. 2.2). Clearly. 



since we do not know iz, we do not know its contour lines dEi either. We do know, 
however, the discrete level characteristics x\ which tell us that, if |Xp ~ Xgl — 1 
for neighbouring sites p and ^, then the contour has to lie somewhere between 
them. In other words, we regard an edge {p, q} of the underlying grid graph as 
intersected by the unknown contour wherever Xp / Xg- Therefore, formula (2.15) 
can be written as 

(2.16) 



where ujpq equals uji from (2.15), if the edge {p, q} is parallel to family Fi. Finally, 



summing (2.16) over all levels / G C and using the coarea formula (2.12) gives an 



estimate for the TV of the continuous signal: 

L-2 

TV(m) « ^ ^c^p^lXp - x[\ = ^^P.kp - u,\ =: TV(u). (2.17) 

The last equality in follows from the fact that X^JXp ~ Xq\ — Ez(Xp ~ Xq)\ — 
\up — Uq\. Of course, the quality of this approximation depends not only on the 
goodness of the Cauchy-Crofton approximation ( |2.15| ), i.e. on neighbour distance d 
and (the size v of) the neighbourhood, but also on the discretization C of R. For a 
relation between the quantized and unquantized discrete TVL^ problems see [15]. 
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2.3.3 Related Approaches 

The authors of [15] define a discrete TV, as any function J : ^ [0, oo] that 
satisfies a discrete coarea formula 

J(u)= / J{x')dl, (2.18) 

and, accordingly, a perimeter for discrete sets A C ^% as Per(A) — J{x^)- This 
definition has two interesting implications. First, such a J has several desirable 
properties, like positive homogeneity, lower semicontinuity and submodularity. 
Second, the straightforward discretization of the TV on a square grid 

/ I V^^l ^ V + - Uj^jY + + l - Uj^jY (2.19) 

^,3 

does not fall into that category, except for the case C — M — {0, 1} where equation 



(2.18) becomes trivial, since x° — u. The TV discretization on the right hand side 



of (2.19) is actually only submodular for £ = B, thus rendering it unsuitable for 
standard graph cut approaches. 

In any case, our derivation of a discrete TV does fulfil ( |2.18| — note that we 



actually already used the discrete coarea formula in (2.17) — , and infinitely many 



more variants can be constructed, namely by using arbitrary neighbourhoods with 



arbitrary positive weights. However, only (2.17) is justified by the Cauchy-Crofton 
formula. 



2.4 A General Version of the Discrete Problem 



After discretizing the fidelity term of (|TyXj in a straightforward way 

\\u- fWl^^n) 



we finally arrive at a discrete version of the TVL"" functional. 

Fa(u) := A^|yp||Mp-/pr +TV(u) 



(2.20) 



(2.21) 



Thus, the combinatorial version of problem (TVL"" ) on an arbitrary lattice reads 

min FJm). (TVr) 

This is the central problem of this paper. As already pointed out in the preceding 
section, functional has the desirable property of being submodular, that is its 
pairwise terms g{up,Uq) = ujpq \up — Uq\ satisfy the inequality 

g{x V 2/) + g{x Ay)< g{x) + g{y), x,y e , (2.22) 

where V, A denote the element- wise maximum and minimum, respectively [41] . 
That (2.22) indeed holds can be shown by using the discrete coarea formula [15]. 

Submodularity implies that multi-label graph cut algorithms can be applied 
directly to solve problem ( TV^"" ) either approximately [TO^ or exactly [34| |2T] . 
However, a discrete TV has (by definition) the even stronger property ( 2.18| , which 
can be used to make solving (TV£^ significantly more efficient by decoupling it 
into binary energies. We refer to ^221 1151 129j on how to do so. 
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2.5 Relations between the Continuous and Discrete Problems 



The continuous TV denoising problem (TVL^) is an infinite- dimensional one. In 



order to be numerically solvable, it must be brought to finite dimensions. This 
is often done by approximating BV{Q) with a sequence of finite-dimensional sub- 
spaces. The general theoretical framework for the discretization of variational 
regularization in a Banach space setting can be found in |48]. In the following 



paragraphs, we briefly examine problem (TV^"" ) within this context. 

A simple subspace of BV{Q) is V(Q%): the class of functions that are piecewise 
constant on the Voronoi cells of lattice Q%. Assume for simplicity u G V{Q%) to 
be quantized, i.e. 



^pX^^ Up^C, (2.23) 



where is the characteristic function of the Voronoi cell of lattice point p. Using 



the coarea formula ( 2.12 ) while observing that level contours of a piecewise constant 



function can be computed by adding side lengths of Voronoi cells, we obtain [TT]: 
TV(^) = E E 2 ^ ^^11^^ "^^1=2 El^^ ^ ^^11^^ - ^^1- (2.24) 

l — O p,q p,q 

The factor | accounts for the fact that every pair of neighbouring cells appears 
twice in the sum. 

Now consider the sampled version of G V{Q%), that is, a discrete image u 
defined on Q%. Let either B — I or B = H, so that the Voronoi cells are regular 
polygons, with side length v. Additionally, let the boundary of Q coincide with 
boundaries of Voronoi cellsj and denote by ^ the nearest-neighbour relation, i.e. 



iff p ^ ^. Then, using (2.24), the discrete functional (2.21) can be 
rewritten as 

Fa{u) = ^^^\Vp\\Up - +CJ^|i^p - Uq\ 

p p^q 

= A||m-/||^ + ^TV(«), (2.25) 

where / is the piecewise constant extension of f = (/p)p, in other words, the 
piecewise constant approximation of data /. Note that for the nearest-neighbour 
relation, the Cauchy-Crofton formula prescribes only one weight oj. We conclude 
that, under the above presumptions, minimization of Fa can be interpreted as 
choosing V{Q%) as approximation space to the continuous problem, that is, solv- 
ing (TVL"") restricted to V{Q%) C BV{Q) only with a different regularization 



parameter A = ^A. Therefore, approximation properties for the space V{Q%) are 
directly passed down to the corresponding discrete problem. 

In [26l [TJ such approximation properties for piecewise constant functions on 
rectangular patches (including V{Qj)) were derived. It was shown that an arbitrary 
u G BV{Q) can be approximated by those simple functions as d ^ 0, only if the 
norm on is used to define the TV (cf. ( |1.1[ )), making it anisotropicjj This can be 
regarded as a decisive factor behind widely observed metrication errors. Roughly 
speaking, the shapes of the constant patches dictate the type of (anisotropic) TV 
which is eventually approximated. This, in turn, causes artefacts in the denoised 
image, since edges that are aligned with grid directions will be "cheaper" than 
others. Similar approximation properties of TV regularization on a hexagonal 
lattice have not been derived yet, but are under investigation in [47] . 



^If this is not the case, the following observation holds up to a negligible constant accounting 
for boundary effects. 

^In [5] these results were extended to general £^ norms, k G [l,co], yet at the expense of 
having to resort either to irregular triangulations of Q or to piecewise linear functions. 
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If larger neighbourhoods are considered, the connection to piecewise constant 
functions does no longer hold. In this case the Cauchy-Crofton weights can be 
expected to reduce the above mentioned discretization effects for any underlying 



lattice (cf. Sec. 2.3.l[ )| | However, they still display one important property even for 
nearest-neighbour grids, which turns out to be useful for experimental issues. An 
easy calculation shows, that the ratio ^ from (2.25) is equal to ^ for both lattices 



A(dl) and K{d' H) and any choices of d! > 0. This implies that, if the discrete 
problem is solved on both grids with the same A, then the regularization parameters 
of the corresponding continuous functionals will be equal as well (A = A^). In 
other words, the Cauchy-Crofton weights ensure that the restored images obtained 
from the two different spatial discretizations correspond to the same degree of 
regularization in the continuous model, and thus can be directly compared (cf. 
Sec. [3}. 

3 Experiments 

Every single one of the following five experiments complies with a fixed procedure, 
which is explained below. 

1. We generate two 8-bit (L = 2^) ground truth images of the same resolution, 
one standard and one hexagonally sampled. For a synthetic test image this 
can be achieved by simply sampling it to two lattices of the same density 



(cf. Sec. 2.1 ). A natural test image has to be resampled to a hexagonal grid. 



since we do not have access to a device producing hexagonally sampled data. 

2. Both ground truth images are contaminated with the same type of noise. 

3. For fixed a (depending on the noise) and a fixed range of reasonable values 
for A (that has been determined before) both resulting versions of problem 



( TYl'^ ) are solved exactly with the sequential algorithm from |22J in combi- 



nation with the max- flow algorithm from [8]. 

4. Finally, steps 2 and 3 are repeated several times to make up for different 
realizations of noise. 

Special emphasis lies on comparing TV regularization with respect to grid topol- 
ogy M% on the one side, to that with respect to Nf and Mf on the other. Apart 
from visually comparing our results, the denoising quality will be measured using 
two different performance figures: first, the distance between ground truth u"^ 
and restored image u divided by the number of lattice points (to make up for the 
fact that the square and hexagonal images will not have exactly the same amount 
of pixels) 

^l|u*-u||i, (3.1) 

which can be interpreted as the mean absolute error per sampling point; and 
second, the ratio of correctly restored pixels 



1 



{pen%:ul = u^}\. (3.2) 



Since the second figure is a less reliable measure of denoising quality, e.g. in cases 
where the restored image suffers from a loss of contrast, it will be only regarded as 
an additional source of information, whereas more emphasis will be placed upon 



(3.1) and visual appearance. At the end of this section the experiments' outcome 



is briefly discussed. 



^Another possibility of reducing metrication errors consists in using higher order pairwise 
terms \up — Uq\^ , (3 > 1 54 . This, however, abandons the TV approach to image restoration 
and is, therefore, out of the scope of this paper. 
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Experiment No. 1 The contrast enhanced Shepp-Logan phantom (Fig. [T]) was 

chosen as test image, and sampled to a square image of size 256 x 256 and to a 
hexagonal one of roughly the same resolution (275 x 238). After adding 60% salt & 
pepper noise, the images were denoised with an data term (a = 1) and usi ng the 
three topologies A/"^, A/"/ and A/"/ and the corresponding weights from Sec. 2.3.1 
The procedure was repeated 50 times for each of several different values of A. This 
setup, a piecewise constant image degraded with impulse noise is perfectly suited 
for TVL^ restoration and can be expected to yield good results. 



While the square grid topologies achieved optimal results in terms of (3.1 ) and 
(3.2 ) for A close to 1, for the hexagonal grid the best choice turned out to be A = 0.9 
(cf. Fig.|7|. This means that "optimally" restored square lattice images are less 
filtered than optimal hexagonal images, since a higher regularization parameter 
means higher data fidelity and less smoothness (cf. the more ragged appearance of 
the square pixel image in Fig. [s] compared to the hexagonal one). If, however, a 
slightly smaller A is chosen for the square image, then certain image details (e.g. the 
tiny ellipses at the bottom) will disappear. Of course, this phenomenon depends 
on the specific test image and, to some extent, the raggedness is also caused by 
the more edgy, i.e. less circular, shape of square pixels. More importantly. Fig. |6] 
indicates that TV regularization on a hexagonal lattice is less prone to metrication 
artefacts. 




Figure 5: Results of Experiment No. 1. Noisy image opposed to denoised images for different 
grid topologies: Affj (centre), Aff (right). The values of A have been chosen to match with 
optimal performance in terms of £^ distance to ground truth (cf. Fig. |7|. 




Figure 6: Results of Experiment No. 1. Details of restored images from Fig, [sj 



Experiment No. 2 Experiment No. 1 was repeated with a different image: 
the cosine from Fig. l2l with a resolution of 270 x 270 pixels. The noise and 
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data term were left unchanged. Since, in this experiment, the differences between 
A/"/ , A/"/ and Mfj were most pronounced (cf. Fig. [9]), the larger neighbourhoods 
J\f}^ and J\f}l were tested as well. Results are presented in Figs, [s] and [9] In 





(bottom right). Values of A have again been chosen to roughly match with optimal 
performance according to Fig. 19] 



contrast to the Shepp-Logan phantom the cosine test image consists entirely of 
smoothly varying intensity changes. TV restoration, however, is known to produce 
cartoon-like images, i.e. minimizers tend to be composed of subregions of more 
or less constant intensity separated by clear edges. This behaviour causes the 
so-called staircasing effect, which is clearly visible in Fig. [s] for both A/"/, where 
the restoration quality is additionally deteriorated by grave metrication errors, 
but also to some extent for Mh- A comparison of N}^ to J\f}^ indicates that, by 
increasing neighbourhood sizes, those effects can be more effectively mitigated on 
the hexagonal lattice. 

It seems appropriate to make some additional remarks on this experiment in 
order to clarify its outcome and especially Fig. [9] The error curves in Fig. [7| and 
even more so in Fig. [9] display a striking feature: at certain distinguished values of 
the regularization parameter A they exhibit significant discontinuities. This pecu- 
liarity of TVL^ regularization, has already been described in T8l|5l]. The authors 
showed that, in general, the data fidelity of minimizers depends discontinuously 
on A, with at most count ably many jumps. This behaviour, which is believed to 
be determined by the scales of image objects that rapidly merge at certain critical 
points, also manifests itself in other quantities, such as our performance figures. 
Fig. 1 1Q| elucidates this phenomenon. 




0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 ' 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 



Figure 9: Results of Experiment No. 2. Averaged distances (left) and averaged ratios of 
correctly restored pixels (right) plotted against A. 




Figure 10: Results of Experiment No. 2. Hexagonally sampled cosine image corrupted with 

60% salt & pepper noise (left), magnified portion of the same image denoised with 
regularization parameter A = 0.906899 (centre), and with A = 0.9069 (right). The image on 
the right, although of the same energy, is significantly less smooth; its TV is increased by 
about 24%, while its data fidelity is accordingly lower. This behaviour corresponds to the 
jump of the black line in Fig. 19] 



Experiment No. 3 The first experiment was repeated with the phantom being 
contaminated by additive Gaussian noise of zero mean and 10% variance. Accord- 
ingly, an data term was employed to remove it. Since the TVL^ model does 
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not preserve contrasts, the number of correctly restored pixel intensities is close to 



zero for reasonable values of A. We therefore only plot the £^ distance in Fig. [12 
This time the differences between grid topologies are less pronounced than in the 
previous experiments. In Fig. ITT] results are presented. 




Figure 11: Results of Experiment No. 3. Noisy image (top left) opposed to denoised images 
with different grid topologies: Affj (top right), Aff (bottom left), Aff (bottom right). The 
regularization parameter A was set to 0.004. 




Figure 12: Experiment No. 3. Averaged distances plotted against A. 
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Experiment No. 4 Experiment No. 1 was repeated once more, but this time 
with a natural image, a 256 x 256 version of the well-known camera man. As stated 
at the beginning of this section, we need to resample natural images to a hexagonal 
grid. A rather simple and straightforward method was adopted in order to do so. 

First, the original image is enlarged by a certain scaling factor c, so that each 
pixel is replaced with a square of c x c pixels (if c is integer). Now, the blown 
up image is covered with hexagonal hyperpixels, from which the pixel values of 
the hexagonal lattice are computed by taking some (weighted) average. If the 
hexagonal image thus generated is supposed to be of roughly the same resolution 
as the original one, the size of one hyperpixel should be close to c^. For this 
experiment we chose c — 7.5 and hyperpixels of size 56 [30], which have also 
been used to create all the hexagonally sampled images in this note. In order 
to avoid smoothing effects as much as possible the median was employed instead 
of the arithmetic mean to compute the new pixel intensities. This results in a 
hexagonally sampled version of the camera man with a resolution of 274 x 240, 



having 0.34% more sampling points than the original image (Fig. 13). 

After adding 60% salt & pepper noise, the images were again denoised with an 
data term. Examples of denoised images are presented in Fig. 14 error plots 
can be found in Fig. |15| The latter figure seems to indicate that the hexagonal 
structure achieved significantly better results for all reasonable values of A. This 
has to be interpreted with care. It is possible, that the better performance figures 
of the hexagonal grid are caused by the grid conversion process simplifying the 
image. 




Figure 13: The two ground truth images used for Experiment No. 4. Original camera man 
image (right) and resampled to the hexagonal grid (left). It should be noted that the camera 
man features many edges that are perfectly aligned to Cartesian coordinates, e.g. contours of 
the camera, tripod and buildings in the background. Naturally those are better captured by a 

square pixel grid. 



Experiment No. 5 Finally, the previous experiment was repeated with another 



natural image and 70% salt & pepper noise instead of 60% (Figs. 16 17) 



Discussion A comparison of the error plots of all five experiments permits two 
remarkable conclusions. First, in every experiment TV denoising on the hexagonal 
grid achieved a smaller minimal distance between restored image and clean 
image. That is, the black curve reaches farthest down in Figs, pf) [9] |12| p^ and|17| 
In the majority of experiments the differences were, however, only marginal. 

Second, if we denote by A"^ the distinguished value, where denoising on the 
hexagonal grid achieves the minimal error, then the hexagonal grid also per- 
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Figure 14: Results of Experiment No. 4. Denoised images with A = 0.9 and grid topologies 

M% (left) a.ndMf (right). 




Figure 15: Results of Experiment No. 4. Averaged distances (left) and averaged ratios of 
correctly restored pixels (right) plotted against A. 



formed better for all (depicted) A < A"^. That is, the black curve lies below the 
others on the left hand side of Figs. [7| [9] [12] [15] and [17] This observation also 
holds true for the ratio of correctly restored pixels (with the only exception be- 
ing Fig. [9|. Note that those A's are more interesting, since larger values result 
in less filtered images, which at least from the viewpoint of visual appearance or 
restoration quality are rather futile (cf. Experiment No. 1 and also Fig. 10). 



4 Conclusion 



The problem of scalar image denoising by means of TV minimization under both 
and data fidelity was considered. We demonstrated how to reasonably 
transform this into a fully discrete (i.e. sampled and quantized) setting on an 
arbitrary lattice. Finally, by conducting a series of five experiments with both 
synthetic and natural images contaminated with different types of noise, we laid 
special emphasis on a comparison of TV regularization with respect to hexagonal 
and square discretizations. 

With the caution in mind that, on the one hand, performance figures can never 
entirely capture the outcome of an experiment and that, on the other hand, eval- 
uation by means of visual appearance is prone to subjectivity and misjudgement, 
it still seems safe to interpret our results as slightly favouring the hexagonal grid. 
This is not only due to the in general better restoration quality in terms of both 
errors and ratios of correctly restored pixels, but also due to its increased capa- 
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Figure 16: Results of Experiment No. 5. Square pixel ground truth image (top left) opposed 
to noisy image (top right) and restored images with respect to topologies Mfj (bottom left) 
and Mj (bottom right). The regularization parameter A was set to 0.9. 




0.6 0.65 

Figure 



17: Results of Experiment No. 5. Averaged distances (left) and averaged ratios of 
correctly restored pixels (right) plotted against A. 



bility of suppressing unwanted metrication artefacts (and, seemingly, staircasing 
effects). Interestingly, the greatest discrepancies between the two discretization 
schemes occurred for the smoothest and most "circular" of test images, the cosine 
in Experiment No. 2, which matches our initial motivation to study hexagonal 
grids. 
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Apart from giving comparative results, our work contributes to supplying the 
(admittedly still modest) demand for image processing tools handling non-standard 
images and, more importantly, makes a case for the significance of spatial dis- 
cretization of continuous problems. Moreover, it suggests new research directions, 
such as generalization of the proposed discretization approach to higher dimensions 
of both image domain and range or an analysis of its approximation properties, 
which is planned to be tackled in future work. 
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